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Abstract. This paper deals with the numerical solution of the Heston partial 
differential equation that plays an important role in financial option pricing, 
Heston (1993, Rev. Finan. Stud. 6). A feature of this time-dependent, two- 
dimensional convection-diffusion-reaction equation is the presence of a mixed 
spatial-derivative term, which stems from the correlation between the two un- 
derlying stochastic processes for the asset price and its variance. 

Semi-discretization of the Heston partial differential equation, using finite 
difference schemes on a non-uniform grid, gives rise to large systems of stiff 
ordinary differential equations. For the effective numerical solution of these sys- 
tems, standard implicit time-stepping methods are often not suitable anymore, 
and tailored time-discretization methods are required. In the present paper, we 
investigate four splitting schemes of the Alternating Direction Implicit (ADI) 
type: the Douglas scheme, the Craig & Sneyd scheme, the Modified Craig & 
Sneyd scheme, and the Hundsdorfer & Verwer scheme - each of which contains 
a free parameter. 

ADI schemes were not originally developed to deal with mixed spatial- 
derivative terms. Accordingly, we first discuss the adaptation of the above 
four ADI schemes to the Heston equation. Subsequently, we present various 
numerical examples with realistic data sets from the literature, where we con- 
sider European call options as well as down-and-out barrier options. Combined 
with ample theoretical stability results for ADI schemes that have recently been 
obtained in In 't Hout & Welfert (2007, Appl. Numer. Math.), we arrive at three 
ADI schemes that all prove to be very effective in the numerical solution of the 
Heston partial differential equation with a mixed derivative term. 
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1 Introduction 



It is a well-known fact that in actual markets the Black-Scholes assumptions are 
violated. The most apparent violation is that the volatility implied from traded 
options, the implied volatility, is not constant but exhibits a strike dependency 
and a term structure. This strike dependency is usually referred to as the im- 
plied volatility smile. The implied volatility smile has serious implications when 
trying to calculate the values of exotic options because they generally do not 
depend monotonically on the volatility. This has inspired numerous authors to 
propose models that, when fitted to market prices, produce stylized effects such 
as the volatility smile. Examples of such models are local volatility processes 
[TTJ HH [57] , jump processes Levy processes jM] and stochastic volatility 

models. A stochastic volatility model that has been particularly successful at 
explaining the implied volatility smile in equity and foreign exchange markets 
is the Heston model. In his seminal paper, Heston [T5j derived an analytic for- 
mula in semi closed-form for the price of a vanilla option, which enables a quick 
and reliable calibration to market prices, especially for liquidly traded vanilla 
options with maturities between 2 months and 2 years [8J. A full and reliable 
calibration of the whole volatility surface can be achieved by adding jumps to 
the asset process [5]. A virtue of the Heston model is that, contrary to e.g. lo- 
cal volatility processes, it displays a realistic behavior in the implied forward 
volatility which is important when pricing forward skew dependent claims |15j . 

The above arguments make the Heston model a prominent candidate for 
valuing and hedging exotic options. Contrary to the Black-Scholes model, to 
date in the Heston model no closed-form analytic formulas have been found for 
exotic option^]. Since no such formulas are available in the literature for any 
but the simplest payoffs, numerical techniques have to be used. 

In the Heston model, values of options are given by a time-dependent partial 
differential equation (PDE) that is supplemented with initial and boundary 
conditions [19] ■ The Heston PDE constitutes a two-dimensional extension to 
the celebrated, one-dimensional, Black-Scholes PDE 6 a . A well-known and 
versatile strategy for the numerical solution of initial-boundary value problems 
for multi-dimensional PDEs is the method-of-lines approach. In this approach, 
the PDE is first discretized in the spatial variables, yielding a large system 
of stiff ordinary differential equations. This, so-called, semi-discrete system is 
subsequently solved by applying a suitable numerical time-stepping method. 

For the semi-discretization of the Heston PDE we consider common finite 
difference formulas. Since the spatial dimension is larger than one, the sizes of 
the obtained semi-discrete systems are in general very large, and standard time- 
stepping methods, such as the popular Crank-Nicolson scheme (trapezoidal 
rule), are often not effective anymore. Accordingly, tailored time-discretization 
schemes are required. 

For the numerical solution of the semi-discrete Heston PDE we shall study 
in this paper splitting schemes of the Alternating Direction Implicit (ADI) type. 

For some recent results see, however, |18| . 
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In the past decades, ADI schemes have been successful already in many appli- 
cation areas. A main and distinctive feature of the Heston PDE, however, is 
the presence of a mixed spatial-derivative term, stemming from the correlation 
between the stochastic processes for the asset price and its variance. It is well- 
known that ADI schemes were not originally developed to deal with such terms. 
In the present paper, we will investigate the adaptation of several important 
ADI schemes to the numerical solution of the Heston PDE with arbitrary cor- 
relation factor p G [—1,1]. The initial and boundary conditions for the Heston 
PDE are determined by the particular option. As test cases, we will consider 
in this paper the pricing of European call options and down-and-out barrier 
options. Through various numerical examples with realistic data sets from the 
literature, combined with ample theoretical stability results that have recently 
been obtained, we arrive at three ADI schemes that all prove to be very effective 
in the numerical solution of the Heston PDE with a mixed derivative term. 

We note that other numerical methods for solving the Heston PDE have 
been considered in the literature. These include schemes based on the finite 
element method gD], the Hopscotch method [H] and fractional step methods 
[131 [25] . Monte Carlo methods for the Heston stochastic process [2] are typically 
employed when the dimension of the problem increases [17 . 

An outline of our paper is as follows. 

In Section [2~T1 we formulate the Heston PDE together with initial and bound- 
ary conditions for European call options. In Section 12.21 we describe the finite 
difference discretization of the Heston PDE. A non-uniform spatial grid is used 
to capture the important region around the strike. In Section l2~3l we study the 
accuracy of the finite difference discretization in various examples of parameter 
sets for the Heston model obtained from the literature. Here the availability 
of Heston's analytic pricing formula for European call options makes an actual 
computation of the global spatial errors possible. In Sect ion HOI we formulate the 
ADI type schemes under consideration in this paper for the semi-discrete Hes- 
ton PDE with a mixed derivative term: the Douglas scheme, the Craig & Sneyd 
scheme, the Modified Craig & Sneyd scheme, and the Hundsdorfer & Verwer 
scheme. Each of these contains a free parameter 9. We discuss the different 
origins of the four schemes and review theoretical stability results that were 
recently obtained in [21[ 122] concerning their application to multi-dimensional 
convection-diffusion equations with mixed derivative terms. In Section 12.51 wc 
perform numerical experiments with all the ADI schemes above, where we ana- 
lyze the behavior of the global temporal errors for each example introduced in 
Section l2~3l As a reference method, we also consider a Runge-Kutta-Chebyshev 
scheme. In Section [2~6l we present numerical experiments for down-and-out call 
options. For these exotic options a closed-form analytic pricing formula has only 
been obtained [57] in the literature if the correlation p = 0. 

Section [3] summarizes our conclusions concerning the four ADI schemes in 
the numerical solution of the Heston PDE with a mixed derivative term. Sub- 
sequently, several issues for future research are discussed. 
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2 Numerical solution of the Heston PDE 



2.1 The Heston PDE 

Let u(s, v, t) denote the price of a European option if at time T—t the underlying 
asset price equals s and its variance equals v, where T is the given maturity 
time of the option. Heston's stochastic volatility model implies [TH1 HZ] that u 
satisfies^ the parabolic PDE 



du < o d 2 u d 2 u , 2 d 2 u du du 



for < t < T, s > 0, v > 0. The Heston PDE dHU) can be viewed as a 
time-dependent convection-diffusion-reaction equation, on an unbounded two- 
dimensional spatial domain. The parameter k > is the mean-reversion rate, 
r\ > is the long-term mean, a > is the volatility-of- variance, p G [—1, 1] is the 
correlation between the two underlying Brownian motions, and r^, Tf denote 
the domestic and foreign interest rates, respectively. In this paper we always 
assume that 2kt] > a 2 , which is known as the Feller condition. 

For a European call option, the payoff yields the initial condition 

u(s, v, 0) = max(0, s - K ) (2.2) 

where K > is the given strike price of the option. Further, a boundary 
condition at s = holds, 

u(0,v,t)=0 (0<t<T). (2.3) 

At the boundary v — no condition is specified. From the assumption 2kt] > a 2 
it follows that this is an outflow boundary. 

As a preliminary step towards the numerical solution of the initial-boundary 
value problem for the Heston PDE, the spatial domain is restricted to a bounded 
set [0, S] x [0, V] with fixed values S, V chosen sufficiently large. The following 
additional conditions at s = S and v — V are imposed for a European call 
option, cf. pT9l [40] : 

^{S,v,t) = e^f* (0<t<T), (2.4) 
u(s, V, t) = se^f* (0<t<T). (2.5) 



Clearly, the boundary conditions (|2.3[) , (|2.5|) are of Dirichlet type, whereas 
(|2.4| is of Neumann type. We take S = 8K and V — 5. This yields a negligible 
modeling error with respect to ()2. !() — ()2.3() on the unbounded domain for a wide 
range of parameter values. 



We assume w.l.o.g. that the market price of volatility risk equals zero. 
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2.2 Space discretization: finite difference schemes 

For the initial-boundary value problem (|2.ip ~ (|2.5p we perform a spatial dis- 
cretization on a Cartesian grid by finite difference (FD) schemes. Here we 
apply non-uniform meshes in both the s- and w-directions such that relatively 
many mesh points lie in the neighborhood of s = K and v — 0, respectively. 
The application of such non-uniform meshes greatly improves the accuracy of 
the FD discretization compared to using uniform meshes. This is related to the 
facts that the initial function (|2.2p possesses a discontinuity in its first derivative 
at s = K and that for v « the Heston PDE is convection-dominated. It is also 
natural to have many grid points near the point (s, v) = (K, 0) as in practice 
this is the region in the (s, u)-domain where one wishes to obtain option prices. 
The type of non-uniform meshes that we employ has recently been considered 
e.g. by Tavella & Randall [37] and Kluge [25]. 

We first define the mesh in the s-direction. Let integer mj > I and constant 
c > 0. Let equidistant points £o < £i < . . . < £ mi be given by 

& = sinh -1 (-K/c) + i ■ A£ (0 < i < mi) 

with 

A£ = — [sinbT 1 ^ - K)/c) - sinh^t-K/c)}. 
mi 

Then a non-uniform mesh = sq < si < . . . < s mi = S is defined through the 
transformation 

Si = K + c smh(£i) (0<£<mi). (2.6) 

This mesh is smooth in the sense that there exist real constants Co, Ci, Ci > 
such that the mesh widths Asi — s,; — s^-i satisfy 

C A£ < As, < Ci A^ and |As l+1 - As 4 | < C 2 (A£) 2 (2.7) 

uniformly in i and mi. The parameter c controls the fraction of mesh points s, 
that lie in the neighborhood of the strike K. In particular, 

Asj w c A^ whenever Sj « if. 

In our numerical experiments we have taken c = K/5. 

We define a non-uniform mesh in the w-direction analogously. Let integer 
77i2 > I and constant d > 0. Consider equidistant points given by rjj = j ■ Ar/ 
for j — 0,1, ... , m 2 with 

At? = — sinh^tV/d). 
m 2 

Then we define a mesh = i'o < vi < . . . < v„ l2 = V through 

Vj = d sinh(^) (0 < j < m 2 ) (2.8) 
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Figure 1: Sample grid denned by J2j|, (J2H) for mi = 30, m 2 = 15, K = 100. 



and write Avj = Vj — Also the mesh (|2.8I) is smooth. The parameter d 

controls the fraction of mesh points Vj that lie near v — 0. It holds that 

Avj ss dAr/ whenever Vj « 0. 

In our experiments we have taken d = V^/500. 

As an illustration, Figure Q] displays the spatial grid defined by (|2.6p , (|2.8D 
for the (small) sample values mi = 30, TO2 = 15 where K — 100. Clearly, there 
are many grid lines near s = K and v = 0. 

We subsequently formulate the FD schemes that we use. Let / : M — > M be 
any given function, let xq < x\ < x 2 < ■ ■ ■ < x m be any given mesh points and 
Axi = Xi — To approximate the first derivative f'(xi), we consider three 

FD schemes: 

f'(xi) w ai- 2 f(xi- 2 ) +Q! i) _i/(^_i) + ai >0 f(xi) , (2.9a) 
/'(or,) « A,-i/(^-i)+A,o/(^)+A,i/(^+i) , (2.9b) 
/'(a;,-) w ji,of(xi) +j itl f(x i+1 ) +ji,2f(xi+2) , (2.9c) 

with coefficients given by 

_ Ai, _ -Axj-i-Azj _ Ai,-i+2Axj 

a ^~ 2 — Ai,_i(Ai,_i+Ai,) ' a »,-l — Aii-iAi; ' a% - a ~ Ai,(Aii_i + Aa:,) ' 

o _ -Azi+i o _ Axj+i-Axj a _ Axt 

Ai.fAi.+Aii+i) ' ~~ AKiAxi + i ' ~~ Ai 1+1 (Ai,+Ai, + i) ' 

-2Aa:^i-Asi^2 _ As;+i+Aa:j+2 — Axs+i 

T*>° — Ai i+ i(Ai i+ i+Ai i+2 ) ' — Ai i+ iAi i+2 ' ~~ Ai i+ 2(Ai, + i+Ai, +2 ) ' 



G 



To approximate the second derivative /"(xj), we deal with the FD scheme 



f"(xi) « 5 h -i fin-i) + 5 h0 f(xi) + 5 iA f(x i+1 ) , (2.10) 

where 

X _ 2 r _ -2 r _ 2 



Next, assume / : K 2 — > M is a function of two variables (x,y). We consider 
approximating the mixed derivative f xy (x,y). Let a^, Ax^ be as above, let 
mesh points yo < J/i < J/2 <•••<?/« in the y-direction be given, and write 
At/j = yj — yj-x- Denote by (3i t u the coefficients analogous to in p. 9b ). but 
then relevant to the y-direction. For the discretization of the mixed derivative 
we use the FD scheme 

d 2 f 1 - 

Q^QZ( Xi ^y^ ~ Pi,kPj,lf{ x i+k>yj+l)- (2-11) 

Clearly, the approximation (|2.11[) can be viewed as obtained by application of 
(|2.9b ) successively in the x- and y-directions. 

The FD schemes above are all well-known. Formulas (12.9b ). (|2. 10|) . (|2. 1 1[) are 
central schemes, whereas (|2.9h ). p. 9b ) are upwind schemes. We note that these 
schemes have previously been applied by Kluge [25j in the numerical solution of 
the Heston PDE. Through Taylor expansion it can be verified that each of the 
formulas (|2.9|) , (|2.10|) , (|2.11[) has a second-order truncation error, provided that 
the function / is sufficiently often continuously differentiable and the meshes in 
the x- and y-directions are smooth, cf. ()2.7|1 . 

The actual FD discretization of the initial-boundary value problem for the 
Heston PDE is performed as follows. 

In view of the Dirichlet boundary conditions (|2.3[) and (|2.5p . the grid in 
[0, S] x [0, V] is given by 

Q = {(si,vj) : 1 < i < mi, < j < m 2 - 1}. 

At this grid, each spatial derivative appearing in (12. lj) is replaced by its cor- 
responding central FD scheme (|2.9b ). ()2.10j) . or (|2.11|) - except in the region 
v > 1 and at the boundaries v — and s = S. 

In the region v > 1 we apply the upwind scheme (|2.9a ) for du/dv whenever 
the flow in the ^-direction is towards v = V. This is done so as to avoid spurious 
oscillations in the FD solution when the volatility-of-variance a is close to zero. 

At the outflow boundary v — the derivative du/dv is approximated using 
the upwind scheme (|2.9b ). All other derivative terms in the ^-direction vanish 
at v = 0, due to the factor v occurring in (|2.ip . and hence, these terms do not 
require further treatment. 

At the boundary s = S the spatial derivatives in the s-direction need to be 
considered. First, the Neumann condition ()2.4|) at s = S implies that the mixed 
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Case 1 


Case 2 


Case 3 


Case 4 


K 


1.5 


3 


0.6067 


2.5 


V 


0.04 


0.12 


0.0707 


0.06 


a 


0.3 


0.04 


0.2928 


0.5 


P 


-0.9 


0.6 


-0.7571 


-0.1 




0.025 


0.01 


0.03 


0.0507 


r f 





0.04 





0.0469 


T 


1 


1 


3 


0.25 


K 


100 


100 


100 


100 



Table 1: Parameters for the Heston model and European call options. 



derivative d 2 u/dsdv vanishes there. Next, the derivative du/ds is directly given 
by (|2.4[) . Finally, we approximate d 2 u/ds 2 at s = S = s mi using the central 
scheme (|2.10p with the virtual point S + As mi , where the value at this point is 
defined by extrapolation using (|2.4p . 

The FD discretization described above of the initial-boundary value problem 
(|2 . 1[) — (|2 . 5[) for the Heston PDE yields an initial value problem for a large system 
of stiff ordinary differential equations (ODEs), 

U'(t) = AU{t) + b{t) (0<t<T), U(Q) = U . (2.12) 

Here A is a given m x m-matrix and b(t) (t > 0), Uq are given m-vectors 
with m = mim2- The vector Uq is directly obtained from the initial condition 
(|2.2p and the vector function b depends on the boundary conditions ()2.3|) - (|2.5p . 
For each given t > 0, the entries of the solution vector U(t) to (|2.12p constitute 
approximations to the exact solution values u(s, v, t) of (|2.ip - (|2.5p at the spatial 
grid points (s, v) £ G, ordered in a convenient way. 

2.3 Spatial discretization error 

In this section we consider four numerical examples and assess the actual con- 
vergence behavior of the FD discretization (|2 . 1 2[) of the Heston PDE defined 
in Sect. 12.21 For any given numbers of mesh points mi, ra-i in the s- and 
w-directions, we define the global spatial discretization error at time t = T by 

e(mi,m 2 ) = max { \u(s l ,v j ,T) - U k (T)\ : \K < s t < < Vj < l} . 

Here u denotes the exact European call option price function, satisfying the 
initial-boundary value problem (|2.ip - (|2.3p for the Heston PDE on the un- 
bounded domain. Next, Uk designates the component of the exact solution 
U to (|2.12p that corresponds to the grid point (sj, Vj). 

Clearly, the global spatial error is defined via a maximum norm. The set of 
asset prices s 6 {\K, and variances v £ (0, 1) in our definition encompasses 
most situations of practical interest. We note that the modeling error, that was 
introduced by restricting the domain of the Heston PDE to a bounded set, is 



8 



Case 1 



Case 2 








Figure 2: European call option price functions u in the four cases given by 
Table □ 



also contained in c(jn\ , vn^) . In our experiments this contribution turns out to 
be negligible. 

For the actual computation of the global spatial errors e(mi, mq) we apply a 
numerical time-stepping schem^l to (|2.12p using a small time step so as to ob- 
tain a sufficiently accurate approximation of U (T) . Subsequently, we employ an 
implementation of Heston's semi-analytical formula [19j to acquire values of it. 
For calculating the single integrals occurring in this formula we use a numerical 
quadrature rule. Numerical difficulties one can encounter in the implementa- 
tion, due to the presence of multi-valued complex functions, have recently been 
discussed in [T]. We adopt the algorithm proposed in loc. cit. where branch 
cuts in the complex plane are correctly taken into account. 

We perform numerical experiments in the four cases of parameter sets given 
by Table Q] Observe that in three of the four cases there is a substantial corre- 
lation factor p. Only in case 4 the correlation factor is relatively small. 

Case 1 has been taken from Albrecher et. al. pTJ, where we have chosen T = 1. 
Case 2 comes from Bloomberg 7 . A special feature of this parameter set is that 
a is close to zero, which implies that the Heston PDE is convection-dominated in 
the w-direction. It is interesting to note that PDEs of this type also arise e.g. in 

3 We used the HV scheme from Sect. E3]with 6 = \ + \ s/Z. 
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Figure 3: Global spatial errors e(2m2,m,2) vs. l/m2 for m-2 = 10,20, . . . , 100 
in the four cases given by Table [TJ 



interest rate modeling, cf. 014]. Values for r^, 77 and T were not specified in |7J 
and have been chosen separately. Case 3 was taken from Schoutens et. al. [35] , 
Here the Feller condition is only just met. Finally, case 4 stems from Winkler 
et. al. 00]. 

Figure [2] displays the exact option price functions u corresponding to the 
four cases of Table[T]on the domain (s,v) £ [0, 200] x [0, 1]. 

Figure [3] subsequently shows for each case from Table Q] the global spatial 
errors e(2m2, ma) vs. l/ma for TO2 = 10, 20, ... , 100. We have taken mi = 2ma 
as it turns out that, for efficiency reasons, one can use much less points in 
the w-direction than in the s-direction. To determine the numerical order of 
convergence p of the spatial discretization, we have fitted in each case a straight 
line to the outcomes for the global spatial errors. Accordingly, in the cases 
1, 2, 3, 4 we found the respective orders p — 1.9, 2.0, 2.1, 2.4. This clearly 
suggests that the FD discretization of the Heston PDE described in Sect. 12.21 is 
convergent of order (at least) two. 

We remark that for the global spatial errors in relative sense, we obtained 
in each case that it is close to 1.0 percent for m 2 = 30 and decreases to ap- 
proximately 0.1 percent for mi = 100. Here we considered asset-variance pairs 
(si,Vj) such that the option value u(si,Vj,T) is always at least 1. 
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2.4 Time discretization: ADI schemes 

Acquiring an effective numerical time-discretization method for the spatially 
discretized Heston problem (|2.12[) is a key step in arriving at a full numerical 
solution scheme for the Heston PDE that is both efficient and robust. 

Let At > be a given time step and let temporal grid points be given by 
t„ = nAt for n = 0, 1, 2, . . .. A well-known method for the numerical solution 
of stiff initial value problems for systems of ODEs 

U'(t) = F{t,U{t)) (0<t<T), U(0) = U (2.13) 

is the Crank-Nicolson scheme or trapezoidal rule. This method defines ap- 
proximations U n to the exact solution values U(t n ) of (|2.13|) successively for 
n = 1,2,3,... by 

U n = U n -i + ±AiF(t„_ x , U n -x) + \AtF{t n , U n ). (2.14) 

In our case of (|2.12[) we have 

F(t, w) = Aw + b(t) for < t < T, w e M m . (2.15) 

Thus, each step (|2. 14|) requires the solution of a system of linear equations 
involving the matrix (/ — ^AtA) where I denotes the m x m identity matrix. 
Since (/ — ^AtA) does not depend on the step index n, one can compute a 
LU factorization of this matrix once, beforehand, and next apply it in all steps 
(|2~T4"]) to obtain U n (n > 1). 

The Crank-Nicolson scheme can be practical when the number of spatial 
grid points m = m\m,2 is moderate. In our application to the two-dimensional 
Heston PDE, however, m usually gets very large and the Crank-Nicolson scheme 
becomes ineffective. The reason for this is that (/ — ^AtA), and hence the ma- 
trices in its LU factorization, possess a bandwidth that is directly proportional 
to minjmi, 7712}. 

For the numerical solution of the semi-discretized Heston problem (|2 . 12[) we 
shall consider in this paper splitting schemes of the ADI type. We decompose 
the matrix A into three submatrices, 

A = A + A x + A 2 . 

We choose the matrix Aq as the part of A that stems from the FD discretization 
of the mixed derivative term in (|2.ip . Next, in line with the classical ADI idea, 
we choose A\ and A2 as the two parts of A that correspond to all spatial 
derivatives in the s- and w-directions, respectively. The r^u term in (|2.1[) is 
distributed evenly over A%, A2. The FD discretization described in Sect. 12.21 
implies that A\ , Ai are essentialljQ tridiagonal and pentadiagonal, respectively. 

Write b(t) from (|2.12p as b(t) = bo(t) + bi(t) + b2(t) where the decomposition 
is analogous to that of A. Next, define functions Fj (j = 0, 1, 2) by 

Fj(t,w) =AjW + bj(t) for < t < T, w e M m . (2.16) 

4 I.e., possibly up to a permutation. 
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Then for F given by (j2~T5)l we have the splitting F = F + Fi + F 2 . Clearly, 
Fq ^ whenever the correlation factor p ^ 0. 

Let 8 be a given real parameter. In the following we formulate four splitting 
schemes for the initial value problem (|2.13|) . Here we assume that F stems 
from a FD discretization of a general 2D convection-diffusion equation with 
a mixed derivative term that is decomposed similarly as above for the semi- 
discrete Heston PDE. All four schemes generate, in a one-step manner, ap- 
proximations U n to the exact solution values U(t n ) of (|2.13|) successively for 
n= 1,2,3,...: 

Douglas (Do) scheme: 

' Y = EZn_i + AiF(t n _i, Z7„_i), 

Y j = Y j _ 1 +9At(F j (t n ,Y j )-F j (t n ^ 1 ,U n . 1 )) (j = l,2), (2.17) 

. U n =Y 2 

Craig & Sneyd ( CS) scheme : 



Y 


= u n . 


_ l + AtF(t n _ 1 ,U n _ 


i). 






Yj 


= Yi- 


r+0A< (Fj (t n , Yj ) 


-Fj{t n . 




(i = 1,2), 


Y 


= Y 


+ ±At(F (t n ,Y 2 )~ 


Fo(t n -i 


,Un-l)), 




Yj 


= %- 




Fj i^n- 




Ci = 1,2), 




= % 











Modified Craig & Sneyd (MCS) scheme: 



Y 


= u n 


_ 1 + AtF(t n _ 1 ,U n _ 1 ), 






Yj 


= Yj- 


i + 9At {F 3 {t ni Y 3 ) ~ F^U-uUn-x)) 


U 


= 1,2), 


% 


= Y 


+ 6 At (F (t n ,Y 2 ) - F (tn_i, CT n _!)) , 








= % 


■f (i - fl)A* (F(t„, y 2 ) - F(t n _i, 


)), 




Yj 




i + #At (F,-(t n , Y,) - Fj(t n _i, U n -i)) 


(i 


= 1,2), 


Un 


= % 









(2.18) 



(2.19) 



12 



Hundsdorfer & Verwer (HV) scheme: 



( Y = E7 n _i + AtFffn-!, C/„_i), 



Yj = Yj-i + 9 At {F 3 {t n , Y 3 ) - Fj(t n -i, U„-i)) 



(J = 1,2), 



< 



Y =Y + 



(F(t n , Y 2 ) — F(t n -i, U n -i)) , 



(2.20) 



Yj = Y^x + 9 At {F 3 {t n ,Y,j) ~ Fj(t n ,Y 2 )) (j 



1,2) 



I U n = Y 2 . 



In the Do scheme (|2.17|) . a forward Euler predictor step is followed by two 
implicit but unidirectional corrector steps, whose purpose is to stabilize the 
predictor step. The CS scheme (|2~T8|) . the MCS scheme (f2~T9|) and the HV 
scheme (|2.20[) can be viewed as different extensions to the Do scheme. They all 
perform a second predictor step, followed by two unidirectional corrector steps. 

Each of the above splitting schemes treats the mixed derivative part F in a 
fully explicit way. Note that in the special case where F — 0, the CS scheme 
reduces to the Do scheme, but the MCS scheme (with 9 ^ |) and the HV 
scheme do not. 

The F\ and F 2 parts are treated implicitly in all four schemes. In every step 
of each scheme, systems of linear equations need to be solved involving the two 
matrices (I — OAtAj) for j = 1,2. Like for the Crank-Nicolson scheme, these 
matrices do not depend on the step index n, and thus one can determine their 
LU factorizations once, beforehand, and next apply them in all time steps to 
compute U n (n > 1). However, the bandwidths of the matrices given by LU 
factorization of (/ — OAtAj) (j = 1,2) are now fixed, i.e., independent of mi 
and rri2- As a consequence, in each of the four splitting schemes, the number 
of floating point operations per time step is directly proportional to m, which 
yields a big reduction compared to the Crank-Nicolson scheme. 

For any parameter 0, the classical order of consistency of the Do scheme - 
for general Fa, F\, F 2 - is just one. An advantage of the CS, MCS and HV 
schemes is that they can attain order of consistency two for general Fq, F\, F 2 . 
Taylor expansion shows that the CS scheme has order two if and only if = h. 
Subsequently, the MCS and HV schemes have order two for any given 0. With 
the latter two schemes, the parameter can thus be chosen to meet additional 
requirements. A virtue of all four splitting schemes is that all internal vectors 
Yj, Yj form consistent approximations to U(t n ). 

The Do scheme can be regarded as a direct generalization of the classical ADI 
schemes for 2D diffusion equations by Douglas & Rachford [12] and Peaceman 
& Rachford [3U] to the situation where a mixed spatial derivative term is present 
in the equation. This generalization was considered by McKee & Mitchell [55] 
for diffusion equations and then by McKee et. al. ^291 for convection-diffusion 
equations. 
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The CS scheme was developed by Craig & Sneyd [TU] with the aim to arrive at 
a stable second-order ADI scheme for diffusion equations with mixed derivative 
terms. 

The MCS scheme has recently been introduced by In 't Hout & Welfert [32] 
to obtain more freedom in the choice of 9 as compared to the second-order CS 
scheme. 

The HV scheme was designed by Hundsdorfer [23] and Verwer et. al. [39] 
for the numerical solution of convection-diffusion-reaction equations arising in 
atmospheric chemistry, cf. also [21] . The application of the HV scheme to equa- 
tions containing mixed derivative terms has recently been studied by In 't Hout 
& Welfert [22 [22]- 

Our formulation of the ADI schemes (|2 . 1 T[) (|2 . 20[) is similar to the type of 
formulation used in [23 . In the literature, these schemes are also sometimes 
referred to as Stabilizing Correction schemes, and are further closely related to 
Approximate Matrix Factorization methods and IMEX methods, cf. e.g. [24] , 

ADI schemes have been mentioned by a number of researchers already for 
the numerical solution of PDE models in finance. Lipton [27J describes an ADI 
scheme for 2D convection-diffusion equations with a mixed derivative term that 
is closely related to the Do scheme with 9 — |. Andreasen [3] H] uses the CS 
scheme with 6 — \ for certain interest rate models without a mixed derivative 
term (when the Do and CS schemes are equivalent) and Randall [32] applies 
this scheme to the Heston PDE with a mixed derivative term. Actual numerical 
results are not presented in these references. The HV scheme has recently been 
considered for the application to PDE models in finance by In 't Hout [20] , where 
an initial experiment in the case of the Heston PDE with a mixed derivative 
term is discussed. As a corollary of our present paper we will find that the choice 
of 9 = 0.3 for the HV scheme [2H] is not optimal with respect to robustness. 

Theoretical stability results for all four ADI schemes - relevant to FD dis- 
cretizations of 2D convection-diffusion equations with a mixed derivative term 
- have been derived in [101 122 1221 I2E1 122] • These results concern uncondi- 
tional stability, i.e., without any restriction on the time step At. The analysis 
in loc. cit. has been performed following the classical von Neumann method 
(Fourier transformation), where the usual assumptions are made that the co- 
efficients are constant, the boundary condition is periodic, the spatial grid is 
uniform, and stability is considered in the (2-norm. 

Currently, the most comprehensive stability results for the Do, CS, MCS 
and HV schemes relevant to multi-dimensional PDEs with mixed derivative 
terms have been obtained in [211 122j . We review the main conclusions from 
loc. cit. relevant to the situation of our paper. 

The Do and CS schemes are both unconditionally stable when applied to 2D 
convection-diffusion equations with a mixed derivative term whenever 6 > i. 
In particular, the second-order CS scheme is unconditionally stable. 

The MCS and HV schemes are unconditionally stable when applied to 2D 
pure diffusion equations with a mixed derivative term whenever 6 > | and 
9 > 1 — |\/2) respectively. Recall that the MCS and HV schemes are of order 
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two for any given 9. At this moment, unconditional stability results for the MCS 
and HV schemes in the general situation of 2D equations, with convection, are 
lacking. It was conjectured [3T] however that the HV scheme is unconditionally 
stable when applied to 2D convection-diffusion equations with a mixed derivative 
term whenever 9 > \ + \ V?>. 

2.5 Temporal discretization error 

In this section we present several numerical experiments for the ADI schemes 
IpTTT]) . ([27T8]) . |[2TT5|) . (|2~20|) which yields important insight into their actual 
stability and convergence behavior in the application to semi-discretized Hcston 
problems (I2.12p with non-zero correlation. 

We define the global temporal discretization error at time t = T = NAt by 

e(N;m 1 ,m 2 ) =max{|[/ fe (T) - U N ,k\ ■ \K < Si < § K, < v 3 < 1} , 

where the index k is such that Uk(T) and Un,u correspond to the spatial grid 
point (si,Vj). Clearly, the global temporal error is defined for the same (s,v)- 
domain as the global spatial error and we also deal again with the maximum 
norm, cf. Sect. 12.31 

Motivated by the theoretical stability and accuracy results discussed in 
Sect. 12.41 we shall consider the Do and CS schemes with 9 = h and the MCS 
scheme with 9 = i. Next, we consider the HV scheme for the two values 
9 = 1- |V2 w 0.293 and 9 = § + ±%/3 « 0.789, to which we refer in the 
following as HV1 and HV2, respectively. We apply all these ADI schemes in 
each of the four cases of parameter sets for European call options in the Heston 
model listed in Table [TJ 

As a reference we also apply the Runge-Kutta-Chebyshev (RKC) scheme. 
This is an explicit second-order Runge-Kutta scheme which has been con- 
structed such that its stability region includes a large interval [— /3,0] along 
the negative real axis. For a complete discussion of this method see e.g. [2~4"Il3"o] . 
The RKC scheme has a free parameter e. Since the Heston PDE contains a 
convective part, we have chosen e = 10, cf. [3S]. In this case (3 ~ 0.34(j/ 2 — 1), 
where v denotes the number of stages of the scheme. Accordingly, in a given 
application to a semi-discretized Heston problem (|2. 12|) with time step At, the 
number of stages is taken as the smallest integer v > + 3 At r[A] where r[A] 
denotes the spectral radius of A. 

Figures [U [5] display for mi = 2m2 = 100 and mi = 2m,2 = 200, respectively, 
the results for the global temporal errors e(7V;mi,TO2) vs. 1/N for a range of 
step numbers N between N — 1 and N = 1000. 

A first observation from Figures [H [5] is that in the cases 1, 3, 4 the RKC 
scheme has global temporal errors that are often comparable to those of the 
MCS scheme. For relatively large At, however, the RKC scheme requires a 
large number of stages v since in all cases r[A] s» 5.1 • 10 4 (if m\ = 2rri2 = 100) 
and r[A] w 2.2 • 10 5 (if mi = 2m2 = 200). As a consequence, in our experiments, 
RKC turns out to be less efficient than MCS for temporal errors larger than 10 -3 
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Case 1 Case 2 

1 




Figure 4: Global temporal errors e(N; 100, 50) vs. X/N in the four cases given 
by Table Q] Schemes: RKC with e = 10 (dotted line), Do with 9 = | (diamond), 
CS with 9= \ (circle), MCS with 9 = \ (grey circle), HV with 6 = 1 - |\/2 
(star) and HV with 9 = | + IV3 (square). 



(if mi = 2m2 = 100) and 10 -4 (if mi = 2m2 = 200). Note that because the 
RKC scheme does produce fair temporal errors for all N in each of the cases 1, 
3, 4, this suggests that the eigenvalues of the corresponding matrices A all lie 
close to the negative real axis in these three cases. In case 2, the RKC scheme 
clearly shows instability for, at least, values TV < 100. We explain this from the 
corresponding matrices A to have eigenvalues in the left half plane that possess 
substantial imaginary parts. This is conceivable, since in case 2 the Heston PDE 
is convection-dominated in the ^-direction, cf. Sect. [231 

For the Do, CS, MCS and HV2 schemes, the Figures 0] [5] clearly reveal 
that in all cases 1-4 the global temporal errors always stay below a moderate 
bound and decrease monotonically with N. This is a very favorable result and 
indicates an unconditionally stable behavior of these ADI schemes in all cases. 
Note that it is a non-trivial result, as it does not directly follow e.g. from the 
von Neumann analysis discussed in Sect. 1^41 
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Case 1 Case 2 

1 




Figure 5: Global temporal errors e(N; 200, 100) vs. 1/JV in the four cases given 
by TableQ] Schemes: RKC with e = 10 (dotted line), Do with 9 = | (diamond), 
CS with 6= \ (circle), MCS with 6 = \ (grey circle), HV with 6 = 1 - |\/2 
(star) and HV with 9 = | + IV3 (square). 



For the HV1 scheme, case 2 shows a peak in the global temporal errors which 
is higher if mi = Ira^ = 200 than if rri\ = 2m2 = 100. We conjecture that the 
HV1 scheme is just conditionally stable, under a CFL condition. 

Subsequent inspection of Figures [H [5] yields for the MCS and HV2 schemes 
in each case 1-4 a convergence behavior of the form C(At) 2 (0 < At < r) with 
constants C, r > that are only weakly dependent on the number of spatial grid 
points m. Hence, the MCS and HV2 schemes show a stiff order of convergence 
equal to two. Clearly this agrees with their orders of consistency, cf. Sect. 12.41 
Remark that this is not obvious, as the order of consistency is a priori only 
relevant to fixed, non-stiff ODEs. 

For the HV1 scheme, we obtain a stiff order of convergence equal to two in 
the cases 1, 3, 4. 

The Do and CS schemes exhibit in all cases an undesirable convergence be- 
havior, with temporal errors that relatively large for modest time steps At. 
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This atypical behavior becomes more pronounced when mi, m^ get larger. 
Though the results are not included in the figures, we remark that the (time- 
consuming) Crank-Nicolson scheme shows a similar behavior. The cause for 
this phenomenon is related to the fact that at s = K the payoff function (|2.2p 
is non-smooth and the Do, CS and Crank-Nicolson schemes do not sufficiently 
damp local (high-frequency) errors incited by this. A remedy for this situation 
is to first apply, at t = 0, two backward Euler steps with step size At/2, and 
then to proceed onwards from t — At with the time-stepping schemes under 
consideration, cf. Rannacher [33J and also e.g. [THl [3T] . By adopting this 
damping procedure, we recover in each case 1-4 a stiff order of convergence 
equal to one for the Do scheme and equal to two for the CS scheme. 

Concerning the implementation, we remark that all codes have been written 
in Matlab version 7.2, where all matrices have been defined as sparse. As an 
indication for the computing times, our implementation of the CS, MCS, HV 
schemes each takes per time step about 0.017 cpu-sec (if mi = 2m,2 = 100) 
and 0.068 cpu-sec (if mi = 2m 2 = 200) on an Intel Duo Core T5500 1.6 GHz 
processor with 1 GB memory. 

2.6 Down-and-out call options 

The FD discretization of the Heston PDE described in Sect. 12.21 is readily 
adapted to more exotic options such as barrier options. We briefly discuss here 
down-and-out call options. If B e (0, K) denotes the down-and-out barrier, 
then the spatial domain becomes [B, S] x [0, V}. Next, the boundary conditions 
(|2.3p . (|2.5|) change to, respectively, 

u(B,v,t)=0 {0<t<T) 

and 

u{s,V,t) = {s - B)e- r f t (0<t<T). 
Finally, mesh points in the s-direction are given by (|2.6[) with 
& = sinh _1 ((B - K)/c) + i ■ A£ (0 < % < mi) 

and 

A£ = — [sinlT^fS - K)/c) - wak~H(B - K)/c)}. 
mi 

The above modifications clearly concern only a few lines in the implementation. 

Lipton [27] derived a semi-analytical formula for the prices of double barrier 
options in the Heston model provided the correlation p = and = Tf. By 
using Lipton's formula with a large value for the upper barrier, the semi-discrete 
Heston PDE (|2.12[) for down-and-out barrier options has been validated in all 
four cases from Table [TJ where we chose the barrier B — 95 < 100 = K and set 
p = 0, Td = rf = 0.03. Moving the upper boundary for s to S = IAK so as to 
reduce the modeling error, again a second-order convergence behavior for the 
global spatial errors is obtained in each case, cf. Sect. 12.31 
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We subsequently applied each of the ADI schemes from Sect. 12.51 together 
with the RKC scheme to the semi-discrete Heston PDE for down-and-out barrier 
options where B, S were taken as above and we considered the four original cases 
from Table [TJ For all schemes it turned out to be advantageous, to a larger or 
smaller extent, to employ the damping procedure at t = to obtain a regular 
behavior of the global temporal errors when the time step At is relatively large. 
The conclusions concerning the observed stability and convergence behavior of 
the schemes are similar to those of Sect. 12.51 In particular, we find that in 
all cases 1-4 the Do scheme has a stiff order of convergence equal to one and 
the CS, MCS, HV2 schemes possess a stiff order of convergence equal to two. 
Next, the global temporal errors for the latter three schemes are always close to 
each other in the cases 2 and 4, but in the cases 1 and 3 the HV2 scheme has 
somewhat larger errors than the CS, MCS schemes. Finally, for the HV1 scheme 
we observe again a peak in the global temporal errors in case 2. In the other 
three cases from Table Q] the HV1 scheme shows a stiff order of convergence 
equal to two. 

3 Conclusions and future research 

Among the schemes discussed in this paper, the MCS scheme with 6 = |, 
used with damping at t = 0, seems preferable for the fast, accurate and robust 
numerical solution of the semi-discrete Heston PDE with arbitrary correlation 
p e [-1,1]. The CS scheme with 9 = \ and the HV scheme with 9 = § + ^\/3, 
both applied with damping, form good alternatives. All three ADI schemes 
show an unconditionally stable behavior combined with a stiff order of conver- 
gence equal to two. For the CS scheme it is essential to apply damping at t = 0, 
whereas for the above HV scheme the error constant can be somewhat larger 
than for the MCS scheme, when damping is used. The RKC scheme as well as 
the HV scheme with 9=1 — lack robustness, as they can have an unstable 
behavior if the volatility-of-variance a is close to zero. Also, the RKC scheme 
appears to be relatively inefficient for moderate error tolerances, which are pre- 
dominant through finance. Finally, the Do scheme shows an unconditionally 
stable behavior but has only a stiff order of convergence equal to one whenever 
the correlation p is non-zero. 

We conclude by mentioning some issues for future research. A main issue 
is a theoretical analysis of the stability and convergence properties of the ADI 
schemes that have been observed in the experiments, cf. Sect. 12.51 Next, a 
study of the performance of ADI schemes relevant to other exotic options in the 
Heston model is of much interest. Finally, ADI schemes can be attractive in the 
numerical solution of other multi-dimensional PDEs from finance with mixed 
derivative terms, e.g. the three-dimensional hybrid Heston-Hull- White model. 
The extension of the ADI schemes (|2.17j) - (|2.20|) to such PDEs is straightfor- 
ward. Positive results on unconditional stability of these schemes in arbitrary 
spatial dimensions, for pure diffusion equations with mixed derivative terms, 
have recently been proved in |22j . 
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